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Abstract 

Wake vortex trajectories and strengths are altered radically by interactions with the ground 
plane. Prediction of vortex strength and location is especially important in the vicinity of 
airports. Simple potential flow methods have been found to yield reasonable estimates of vortex 
descent rates in an otherwise quiescent ambient background, but those techniques cannot be 
adjusted for more realistic ambient conditions and they fail to provide satisfactory estimates of 
ground-coupled behavior. The authors have been involved in a systematic study concerned with 
including viscous effects in a wake-vortex system which is near the ground plane. The study has 
employed numerical solutions to the Navier-Stokes equations, as well as perturbation techniques 
to study ground coupling with a descending vortex pair. 

Results of a two-dimensional, unsteady numerical-theoretical study are presented in this 
paper. A time-based perturbation procedure has been developed which permits the use of 
analytical solutions to an inner and outer flow domain for the initial flow field. Predictions have 
been compared with previously reported laminar experimental results. In addition, the influence 
of stratification and turbulence on vortex behavior near the ground plane has been studied. 
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Introduction 

It is known that aircraft trailing vortex wakes can cause serious loss of control when following 
aircraft encounter them. The hazard is more severe near the ground because of the limited time 
and space available to maneuver the aircraft. Therefore, prediction of wake vortex trajectories 
and strengths is especially important for effective airport flight control. The work reported here 
has been focused on the viscous interaction between vortex wakes and the ground plane. 

Potential theory representation of a pair of counter-rotating vortex filaments above an 
infinite plane yields reasonable estimates of wake vortex descent rates in an otherwise quiescent 
atmosphere, when the wake is sufficiently far from the ground (Saffman 1979). By assuming the 
vortices could be treated as small core sized vorticity spots above a very thin ground boundary- 
layer, Liu and Ting (1987) determined that away from the small cores and the thin boundary-layer 
region, the flowfield obeyed the Euler equations. From that perspective, it is very inefficient and 
expensive to study the trailing vortex problem using Navier-Stokes solvers for the entire flow 
field. However, the methods of Liu and Ting (1987) are only valid when the vortical spot is 
sufficiently far away from the ground relative to the core size, so that vortex-core interactions 
with the ground-plane boundary layer can be neglected. Hence, Euler equations can be used to 

model large portions of the computational domain, but that region must be adjusted continuously 
near the ground plane. 

It should be stated that without ground effect, even in more realistic atmospheric conditions, 
an approximate model developed by Greene (1986) demonstrated surprisingly good agreement 
with experimental data. That method was based upon empirical engineering approximations. 
But from the authors’ experience, the extension required to include ground effects was not 
compatible with Greene’s model. 
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Experiments to establish the features of the flowfield induced by a single vortex near the 
ground were carried out by Harvey and Perry (1971) in a low-speed wind tunnel. The primary 
objective of those experiments was to explain the cause of the vortex rebound phenomenon 
which is considered to be a significant feature of the ground effect and which can influence 
terminal flight conditions. They inferred that rebound was caused by separation of the ground 
boundary-layer flow underneath the vortex. They argued that boundary-layer separation produces 
a secondary vortex whose development makes the primary vortex rise. Later Barker and Crow 
(1977) observed rebound for a vortex pair, generated in water, approaching either a free upper 
surface or a rigid horizontal plane immersed in the water. They asserted that the phenomenon 
could be attributed to the effect of finite vortex core size. Recently, experiments were designed 
by Liu and Smsky (1990) to minimize the sidewall effects that can mask the ground effects. 
From their dye visualization results in water, they identified the emergence of secondary, counter- 
rotating vortices outboard from the vortex wake near ground level. They determined that as soon 
secondary vortex began to form, rebound of the main vortex was initiated. Specifically, 
the primary and secondary vortices form a vortex pair that moves upward. Those results have 

confirmed essentially the scenario of secondary vortex generation suggested by Harvey and 
Perry (1971). 

Many theoretical studies of a vortex pair in ground effect have been reported during the last 
two decades. Research reported by Bilanin, Teske and Hirsh (1978), Saffman (1979) and Peace 
and Riley (1983) disagreed with the finite core size hypothesis of Barker and Crow (1977). 
Saffman (1979) showed, within the framework of inviscid theory, that it was not possible to 
explain the rebound phenomenon by finite core size and that the wallward velocity component 
cannot change sign. Navier-Stokes computational results for a trailing vortex pair near the ground 
were performed by Bilanin, Teske and Hirsh (1978) for both laminar and turbulent conditions. 
They demonstrated that the rebound did not occur unless the viscous, no-slip boundary condition 
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was applied. One of the important conclusions drawn by them was that the proximity of a ground 
plane reduces the vortex hazard by scrubbing. That is, the vortex pair separates or spreads and 
interacts viscously with the ground thereby reducing its strength more rapidly. A numerical study 
of the effects of stratification and wind shear on the evolution of aircraft wake vortices near the 
ground was also included in the report by Delisi, Robins and Fraser (1987), who found that 
both effects reduce the extent of vortex rebound. As they stated, turbulence was not included 
formally in their computational model. But some ad hoc small scale damping was added to the 
equations included in the numerical model. They thought this simple “turbulence model” was 
adequate to predict vortex migration. In addition, a mixed no-slip/slip boundary condition was 
invoked on the ground plane which needed an empirical adjustment for different flow cases. 

In the work which follows, a two-dimensional, unsteady numerical- theoretical study will be 
reported. A time-dependent double series, asymptotic expansion in terms of Reynolds number 
and time has been used for the initial flow field. The computational domain and grid spacing 
have been chosen carefully to avoid using approximate boundary conditions and to get proper 
resolution. Subsequently, the influences of stratification and turbulence on vortex behavior near 
the ground plane have been shown. 

Computational Approach 

An unsteady vorticity, stream-function formulation has been used in the numerical calcula- 
tions. Invoking symmetry of the vortex pair permits the computations to be restricted to the first 
quadrant (0 < x < oo, 0 < y < oo). The system of equations were made dimensionless using 
the initial vortex half-span, s 0 , as the characteristic length and the initial circulation, To, along 
with the fluid density, p 0 , to formulate the other dimensionless variables. The characteristic flow 
Reynolds number is Re = Tq/uq, where i/ 0 is the kinematic viscosity. Characteristic velocity, 
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and pressure « given by r„/s„, ,j/r 0 and ^,4. respective*. Since .here is on, y one 
component of vordcity. f . the governing vorticity M d stream-funenon (« equauons are: 




(i) 


and 


V 2 


0 = -C 


respectively. Here 


( 2 ) 


u = ui + vj 


( 3 ) 


and 


dv du 

= ai~d^’ with 


u 


The coordinate system is shown in Figure 1. 


drp_ 

dy' 


and v = — 


di> 

Ih ‘ 


( 4 ) 



Figure I. Coordinate system employed in this study 



Since the flow is for the most part inviscid, with an unsteady, viscous ground-plane boundary 
layer and a small viscous vortex core region, problems were encountered in starting the numerical 
calculations. The initial velocity field problem has been addressed previously by Peace and 
Riley (1983), but they used a single series expansion in time, scaled by Reynolds number 
(i/R e ), to start their asymptotic solution. Limitations of the single series expansion for unsteady 
boundary-layer flows have been discussed by Nam (1990). Considerable effort was devoted to 
developing an appropriate initial velocity field which accommodated the ground effect region 
in the numerical simulation without producing non-physical, numerical start-up transients. A 
forthcoming contractor report (Ash and Zheng 1991) contains a more detailed derivation of 
the asymptotic expansion formulations, but the essential features of the approach are developed 

below. In addition, the computational domain transformation procedure and boundary condition 
specifications are discussed. 

Flow Field Initialization 

Since a goal of this study was to extend the Reynolds number range over which vortex- 
ground plane interactions could be modeled, viscous effects were anticipated in both the vortex 
core and boundary layer regions long before any interactions between the vortex core(s) and 
boundary layer occurred. Furthermore, a fine numerical grid was required for both the vortex 
region, including its path of descent, and the ground-plane boundary layer. Both requirements 
demanded an initial velocity field which was free of anomalous velocity gradients. 

The Oseen (1911) vortex is an exact solution to the Navier-Stokes equations for the diffusion 
of a vortex filament into a viscous region of infinite extent. For any time greater than zero, the 
Oseen vortex includes viscosity, while at t = 0, it is an inviscid point vortex. Hence, placing a 
pair of those vortices at ±x 0 , yo at t = 0, is equivalent to placing a pair of potential vortices at 
those locations, but then allowing viscous effects to occur immediately after placement. Oseen 
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vortex solutions do not include either non-linear coupling of the vortex pair or the viscous 
influence of a ground plane. Adding a mirror image pair of Oseen vortices (at ±x 0 , - yo ) 
can be used to initiate a ground plane interaction. Since those vortices proceed immediately to 
viscous flows, they appear to be a more realistic starting flow than that of a pair of potential 
vortices in any practical numerical grid. Hence, analytic perturbation methods can be used 
to predict the flow during the initial time interval when local velocity gradients create severe 
problems for numerical techniques. 

In the ground plane boundary-layer region, the initial stream-function ip(x,y, t ) is assumed 
to take the form: 


= 2e[^i(x, T),t) + trj> 2 (x , T) , t) + ■ ■ ■] ( 5 ) 

where e = JtJWt and 77 = y/2e, with 

OO 

rl> l {x,r 1 ,t) = Y J tP <t> i p {x,T),t) ( 6 ) 

P = 0 

This stream- function must satisfy the no-slip boundary conditions. 

In the outer flow, it is assumed that the stream -function can be represented initially by 

ty(x,y,t) = (x,y,t) + + ... ( 7 ) 

This stream-function must satisfy the initial potential vortex requirements of both a symmetry 
plane and a ground plane. Application of the two-dimensional, incompressible Navier-Stokes 
equations to the two stream-function series and matching the expansions using van Dyke (1976) 
type matching procedures yields 

= 2e(V>i + exp 2 ) + 0[e 3 ] (8) 

where 

</>. = Ob/01 + t(j£fn + O^J„) + 0 [t 2 ] 


( 9 ) 
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/l2 = 


5 12 


^2 = f/l501 + + ^0^512 + +0[< 2 ] 


( 10 ) 


/ 01 = rj erf(rj) + -L e ^ - -ij 


(ID 


37T 1 / 2 


fn-^rxll- e-" 2 ) - -t? erfc(r}) + ^9^ 


( 12 ) 


6^1/2 ( X e " er / c (^)) + i7^T7I er /(^) - ^( 2 + ^)v 3 erfc(ri) 


~ ^^Mv) + ^Ve 2v 2 - ^ trfc\rj) + ^ 

F (2 “ 3^) ,?er/c(,?) " dW r/(7?) + (Jtj + 


2tt1/2 9^3/2 


1 — e~ v 


(13) 


5oi - 9 + -T: 1 ^r) 2 erfc(T)) - l-rje 1,2 - -x 1 / 2 


erf( 9 ) 


(14) 


5ll = 


ir 1 / 2 


(5 3 + 9*)erfc(r]) + ^ J (-J? 2 + V^^drj 


(15) 






7T 


1/2 


(7 + ^ 3 ) er /(^ + ^( l - f 7 2 ) e -’' 2 


+ 7r ^ ( 2^ 7 ? 3 ^( er /(5)) 2 + - 2rj 2 ^e ,) erf(rf)-^Tje 

- tt^t, erf(rj) - ~erf{ V ) + ? - ie"’ 2 


(16) 


5!3 - 8 tt 1/2 (V + 3 5 4 j er/cfa) + ^J q (V + e'^dT) 


(17) 
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It is noted that Uq(x, t) and U\{x,t) are the outer flow representations of the vortices along 
the ground plane (y = 0). These solutions are restricted to small times and are given by 


Uo{x, t) = - 


Vox 


7T 


(x 1) + j/q (^ + 1) 2 +j/q 


(18) 


and 


V i(x 


2 r 


dU o 
dx 


dj 

X ~ £ 


(19) 


which are the same forms used by Peace and Riley (1983). 

In effect, four Oseen vortex solutions, with Cartesian velocity components given by 

Uc{x,y, t) = (y - y 0 )[$(x,y, t\ — 1 , j/o ) - $(x, jM; l,j/o)] 

+ (y + yo)[$(x,y,t\ 1 , -J/o) - $(x,y,t; — 1 , — J/o)J 


( 20 ) 


and 


Vc(x, y, t) = (x- l)[$(x, y , t ; 1, y 0 ) - $(x, y, t ; 1, -y 0 )] 

+ (x + l)[$(x,j/,f;-l,- yo ) - $(x, y, f; — 1, r/ 0 )] (21) 


where 


1 1 - e -*«[(*-«) 2 +(jf-0) 2 ]/4< 


( 22 ) 


2?r (x - a) 2 + (y- P) 2 
have been employed in the outer flow, in the vicinity of x = 1, y = y 0 . Again, these functions 

are restricted to small times, where they do not alter the asymptotic matching conditions. 


To summarize, the asymptotically generated velocity field, which constitutes the initial 
velocity distribution in the numerical scheme, is a combination of solutions (7), (8), (20) and (21). 
Small values of time, compared to the Reynolds number, have been employed. The asymptotic 
solution time level selected to generate the initial velocity field depended on the circulation 
Reynolds number, but it was the maximum time allowable by the finite term approximations 
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to the infinite series expansions. The asymptotic solutions were thus used to produce an initial 
velocity field whose local gradients could be handled by the numerical grid, but whose series 
representations did not represent a significant amount of computational overhead. 

Domain Transformation and Boundary Conditions 

A moving grid was considered for this study but was discarded subsequently, due to the 
complexity of the vortex trajectory in ground effect and to uncertainty in the viscous zones 
requiring fine grid resolution. It was deemed more reasonable to pack grid points adjacent to 
the symmetry- and ground-planes and allow the vortex system to move with respect to the grid. 
This was possible because the nominal regions where viscous effects occurred were known to a 
first order approximation. It is noted further that the elliptic character of the incompressible flow 
field mandates enforcing the boundary conditions at the infinite limits of x and y or alternatively 
developing rigorous boundary condition approximations within a finite domain. 

Bilanin et al (1977) and Ting (1983) studied the far field boundary condition problem. They 
used the far field expansions of Poisson integrals which could be employed in a finite subdomain 
of an unbounded fluid. The exponential decay laws for voracity distributions were required for 
the convergence of the far field expansions (Ting 1983). The expansions were developed further 
by Ting (1983) using integral invariants. Unfortunately, neither the decay laws nor the integral 
invariants exist when a no-slip boundary is present because the vorticity disaibution in the viscous 
ground-plane boundary layer cannot be estimated a priori. Here a coordinate mapping has been 
used, employing a simple exponential transformation which yields densely packed grid points 
near the ground plane and stretches the grid as the infinite limits are approached. The mapping 

X=a(l-e~ bz ) 

Y = c(l - e- 4 *) 


( 23 ) 
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transforms 0<x<oo, 0 < y < oc into the finite domain 0 < X < a, 0<K<c. The 
corresponding spacing intervals are 


A X % abe~ bx Ax 

AY ^ cde~ dy Ay (24 ) 

Hence the boundary-layer and vortex core resolution can be controlled by adjusting the arbitrary 
constants a, b, c and d. 

Since the ground plane spacing is compressed automatically via the coordinate mapping, 
resolution requirements are more severe in the far field. In addition, vortex core resolution in 
the vertical direction is most sensitive at start-up (y = yo ), while horizontal resolution is least 
accurate at the end of the numerical simulation (when the vortex core is at the greatest horizontal 
distance from the symmetry plane) - say x,. Since the vortex dilates as time increases, selecting 
an appropriate grid spacing initially produces a conservative spacing level at later times. 

As an example, if the initial vortex core radius is r c = 0.2, then we would require that 
Ai/max = 0.02 (= r e /10) at y 0 , and that Ax max = 0.1 (= r c /2) at x,. The Ax max spacing will 

be more compact relative to the vortex core due to dilation beginning from t = 0. However, 
the spacing interval is given by: 


AX < 0.1 abe~ bxi 


and 


AY < 0.01 cde~ ay ° (2 5) 

Since the number of grid increments in the X-direction must satisfy MAX = a, while 

NAY = c. 


M > 


10e 6 *' 


b 
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and 


> SO^* 

- d 

(26) 

The minimum number of grid points for acceptable resolution 
d = l/y 0 and is given by 

occurs when 6 = l/x, and 

M > 10x/e 


and 


N > 50y o e 

(27) 

The boundary conditions in the transformed domain are: 


C(o,r,o = o, 0(0, Y,t) = o, 

(28) 

II 

o 

II 

© 

(29) 

«^o,«) = - c V 0 (A-,o,o, 

(30) 

■ip(X,0,t) = 0, 

(31) 

C(X,c,t) = 0, and xp(X,c,t) = 0. 

(32) 

The discretized expression for the ground-plane vorticity, Eqn.(30), 
(1972), and is given by 

was developed by Roache 

C(*. O,t)^-20 n c 2 d 2 /AK 2 

(33) 
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where the subscript, n, indicates the grid locations adjacent to the ground plane. 

In this study, the vertical grid spacing was fine enough to permit implementation of the 
no-slip boundary condition on the ground plane without exaggerating numerically the viscous 
interaction encountered previously by Delisi, Robins and Fraser (1987). 

While the velocity boundary conditions appear to be over-specified (Anderson 1986), they 
are compatible with the velocity components in the transformed domain, given by 

dxb 

u = ^.d.(c-Y) 
dx b 

v = ~dX’ b ( a - x ) . (34) 

which can be easily verified. These boundary conditions also satisfy the integral conditions 
developed by Anderson (1986). 

An alternating-direction implicit (ADI) scheme was used to solve the vorticity transport 
equation (1), with the upwind flux-splitting method applied to the convection terms and central- 
differencing for diffusion terms. An efficient Poisson solver (Swarztrauber and Sweet, 1979) 
was used ,o solve Eqn.(2). Viscous flow in a driven cavity was chosen as a test problem for 
that computational scheme to explore the capability of capturing the secondary-vortex evolution 
phenomena at high Reynolds number. The systematic numerical studies of this problem are 
well-documented in Bozeman and Dalton (1973) and Rubin and Harris (1975). 

Stratification Effects 

The present analysis has employed the Boussinesq approximations to model buoyancy effects. 
Specifically, the local density has been assumed related linearly to temperature via 


P = po[l - a(T - To)] 


(35) 
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where a is the volumetric coefficient of thermal expansion and a<l. Hence, the density 
departure from the reference value can be represented in dimensionless form as 

p = (p - Po)/po = -a(T - To) ( 36 ) 


The complete Boussinesq model neglects all density variation effects except the body force 

term in the momentum equation so that the conservation of momentum equations result in the 
modified vorticity transport equation: 


where 



1 dp 1 

Tl~di + Te 


v 2 c 


(37) 



(38) 


For the case of a linearly stratified ambient fluid, the dimensionless temperature, 9, defined 
by 

8 = (T - Tq)/Tq ( 39 ) 

is assumed given in the quiescent, ambient state as 


0 = fiy 

Hence, the ambient density departure is given by 

p = -a(3yT 0 

or 


(40) 


(41) 


9 = -p/aTo 


(42) 



Similarly, the density and temperature fluctuations are related such that 
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and 


P = P- P 
O' = 9-9 


(43) 


P = -aT 0 0' 


(44) 


Neglecting viscous dissipation, the appropriate conservation of energy equation, invoking 
the Boussinesq approximation is: 

80 ( i 

aI + { uV ) e = ^p r ^ (45) 

which can be written in terms of the density fluctuations (through Eqn. 44) for the stably stratified 
case as 


dp' . dp' dp' dp 1 

~jrr + — t- v— — f- v— 

ot ox dy dy RePr 


vy 


(46) 


Restricting attention to cases where RePr > 1, we can neglect the diffusion terms and the 


conservation of energy equation becomes: 



(47) 


Alternatively, buoyancy effects could have been introduced through the compressible conser- 
vation of mass equation. By expanding the conservation of mass equation in terms of the volu- 
metric coefficient of thermal expansion, a, the incompressible continuity equation and Eqn. (46) 
evolve from the zeroth and first order groupings in the a expansion. Hence, the numerical 
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simulation is modified via Eqns. (37), (41) and (47) when stratification effects are included. The 
computational scheme was easily modified to include stratification since the same grid was used 

and the density fluctuation equation could be implemented using the same procedures employed 
for the vorticity transport equation. 

It is noted that the dimensionless parameter, F v , which appears in the vorticity transport 

equation (37), is related to the Brunt- Vaisala frequency. That is, if the dimensionless density 
gradient is defined by 



then the Brunt- VdisSlS frequency, N, is given by 


N 2 _ _ 9 _dp = T 2 n 2 
•So dy 


(48) 


(49) 


Stratification effects represent additional difficulties because of the characteristic time in- 
troduced via Brunt- Vaisala or density induced oscillations. If s 2 Q N/r 0 (= -fc) > 1, buoyancy 
effects are significant and the characteristic time Cs 2 /r 0 ) is large compared with the time in- 
terval over which density induced oscillations occur. The complication arising from these two 
characteristic times has been examined in detail by Hirsh (1985). For more realistic physical 
problems, F v is typically quite large and the stratified density gradient (through n) is not large 

enough to result in large values of sIN/Tq. Consequently, a limiting test case is when the two 
time scales are equivalent (^Ayro % 1). 

Simulations were run with aJiV/r, = 1. but the density effects were so large that major 

vortices of opposite sign were generated very rapidly and the flow quickly became unstable 
numerically. 
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Ambient Turbulence 

Inclusion of turbulent effects is an important element of this study. However, due to the 
approximate nature of turbulence models, initial work has restricted attention to limited types of 
turbulent processes and simple models. I, is assumed here that atmospheric turbulence is being 
generated by axial (z-direction) wind shear. Since stratification effects are under investigation, 
it should be noted that buoyancy induced turbulence effects are usually small in the immediate 
vicinity of the ground, when compared to wind-driven fluctuations and have been neglected. 
However, some very important types of buoyancy-driven turbulence phenomena are known to 
occur with potentially serious consequences (e.g. microbursts), but those processes are too 
complicated to be modelled reliably at this time. These catastrophic-type, bouyancy-driven 
turbulence effects supercede wake-vortex prediedon requirements and should be investigated 
separately. In addition, turbulence is generated by the wake vonex structures themselves. 

Those turbulence generation processes are also outside of the turbulence modelling capabilities 
employed in the present study. 

The dimensionless vorticity transport equation can be written in index notation as 

DQ, dui 1 „ 

Dt ~ a ’dT t = n ' (50) 

By assuming that the vorticity components can be decomposed into (slowly varying) time mean 
and (rapidly varying) unsteady parts, i.e. 

= n, + ^2) 

with 

U, = Uj -f- uj- 



then the Reynolds averaged vorticity transport equation can be written: 
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dQj 

dt 




dxj 

_d_ 

dx, 




dxj 



( 52 ) 


Here, we restrict attention to mean flows which are steady and parallel in the z-direction, 
i.e. w = w(y) so that 


u = u(x,y,t) ( 53 ) 

and 

v = v(x,y,t) 

These are essentially two-dimensional flows, but permit the existence of two mean vorticity 
components. That is. 


and 


= dw 

n, = d^ = n,M 


= ((x,y,t), 


( 54 ) 


with 


f^2 = 0 


Consequently, the governing equation on ( becomes: 



( 55 ) 


where the u'ft' terms are the terms which must be modelled. 
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We have already stated that atmospheric turbulence has not been coupled with vortex- 
generated turbulence in this study. Essentially, it is assumed that the aircraft vorticity field 

behaves like a passive scalar, immersed in atmospheric turbulence. In that context, we have 
assumed that 


tt'u' = -ci 9 A 


dxj 


where q 2 is the local turbulent kinetic energy, 


(56) 


while A is a turbulent length scale and c, is the turbulence modelling constant. It is assumed 
that q and A are not altered by vortex interactions. 

In order to utilize Eqn. (56), it is necessary to model q and A. Bilanin, Teske and Hirsh 

(1978) have employed a second order closure model, using similar nomenclature, to study vortex 

wake decay. When the wake vortex velocity field did not contribute to the turbulence, it was 

possible to show that q was constant both near the ground plane and in the far field. In addition, 

the von Karman constant, «, can be used to model the characteristic turbulent length scale near 
the ground as 


Since A should be constant away from the ground, it was necessary to model the intermediate 
zone. Here, we have assumed that 


A«c 2 (l-e-»/«) (59) 

which matches Equ. (58) for small y. von Karman's constant has been taken as 0.4 in this 
s tudy, leaving C 2 as the remaining arbitrary constant. 



Employing Eqn. (56) in Eqn. (55) yields: 
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dy 

+c\q- 


dx 2 


5 (k^\ 

%Vai) <«> 

where A is given by Eqns. (58) and (60) and c\, c 2 and q are evaluated using experimental data. 


Results 


Based upon previous numerical studies, it was determined that the vortex core centers could 
be placed at x 0 = 1, y 0 = 2 (and x 0 = -1, y 0 = 2, from symmetry) to start the simulations. That 
vertical distance (y 0 = 2) was deemed close enough to the ground plane to produce detectable 
coupling effects after moderate time intervals after start-up, but it was also far enough from 
the ground plane to enable the vortex flow field to establish itself prior to strong ground-plane 
interactions. The initial vortex core was assumed to have a core radius, r c , of 0.2 and a 150x300 
grid was employed in the numerical simulations reported herein. Calculations were performed 
using a Cray O computer at NASA Langley Research Center. 

Limited experimental data were available for numerical validation studies. Only the experi- 
ments of Liu and Smsky (1990) were used. Their vortex flows were produced using an NACA 
0012 wing model and the estimated circulation based Reynolds number (T 0 /y) was 7,650. Their 
experiments did not investigate stratification effects near a ground plane. 

For the unstratified, validation studies, circulation Reynolds numbers of 1000, 7,650 and 
75,000 were simulated Figure 2(a) represents the trajectories (*(<), y (t)) followed by the 
three simulated vortices, along with the measured trajectory of Liu and Smsky (1990). The time 
histones of x(<) and y(t) for the four vortex cases are shown in Figures 2(b) and (c), respectively. 
The agreement between the numerical simulation and the experiment is quite good. 
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In order to assess the influence of Reynolds number, stratification and turbulence on the 
vortex hazard, some measure of hazard strength was required. However, since the computational 
domain is an unbounded quadrant, overall or global measures of circulation or velocity levels 
appeared to be of little value. It was finally decided that circulation and kinetic energy histories in 
the aircraft approach zone would be meaningful. The zone was selected somewhat arbitrarily to 
be the area bounded by -2 < x < 2 and 0 < y < 3, which would span a typical runway entrance. 
The velocity components and vorticity were compared at each grid point within the right half 
of that area(0 < i < 2,0 < y < 3) at each time level. Subsequently, zonal circulation, T(2x3), 
was calculated by integrating the vorticity over the half area and the instantaneous kinetic energy 
within the zone, E(2x3), was compared by a similar integration of ^1. Those histories are 
shown in Figure 3 for the extreme Reynolds number cases(1000 and 75,000). The experiments 
of Liu and Smsky did not report circulation and energy history data of this type. 

Data were not available for comparison between simulated vortex flows in a stably stratified 
ambient environment with ground effect. Stratification effects were rested for a circulation-based 
Reynolds number of 1000, at dimensionless Brunt- VaisStt frequencies (n /F.) of N = 0.05 and N 

= 01 - Th ' prediC,ed rcsults for tatar flow case are displayed in Figure 4, along with the 
unstratified case (N = 0), for comparison. The trajectories, lateral and vertical position histories, 
along with circulation strengfl, r ( 2x3) and kinetic energy E(2x3) histories are displayed in 
Figures 4(a) through (e), respectively. 

The background turbulence model, discussed in the previous section, was incorporated in the 

simulations for a circulation-based Reynolds number of 75,000. A global turbulence constant, 

C, defined by 

C = c\C2q 



22 


is all that is required to model this eddy-viscosity type turbulence (since c 2 was chosen as 0.5) 

and values of C = 0.001 and 0.1 were employed in the simulations. The combined effects of 

Reynolds number, turbulence and ground effect on vortex motion and strength are shown in 
Figures 5. 

Discussion 

While the laminar flow simulations (Re = 1000) are less realistic, in terms of aircraft vortices, 

they are less ambiguous, in terms of effects of turbulence models and numerical uncertainties. 

Consequently, vortex rebound and stratification effects, derived from fundamental phenomena, 

can be discussed with more certainty for those flow cases. To that end, stratification effects 

in the vicinity of the ground plane have altered vortex trajectories rather remarkably, as shown 
in Figure 4 . 

Refening to Figure 4, it can be seen that the vortex appears to literally try to “fall back 
down hill” in both stratified cases, even though the initial descent and rebound trajectories nearly 
coincide with the unstratified case. The mechanism which is responsible for that effect can best 

be explained by comparing the computational flow visualization results from the unstratifted 
case (N = 0) with the stratified case (N = 0.1). 

Figures 6 and 7 are color panels comparing the vorticity distributions at four time levels 

for an unstratified flow (Figure 6) and a stratified flow (Figure 7). Streamfunction contours are 
compared in Figures 8 and 9. 

The stratification effect on vortex trajectory is most easily understood by looking at the 
density distribution. Figure 10, shows density departure contours (from Eqn. 41, via p\ given 
in Eqn. 43) a. four different times. There, it can be seen that relatively higher density fluid 
is pulled from the ground plane around the primary vortex, where it tends simultaneously to 
compress the vortex and cause more rapid vortex deceleration due to increased inertia. While 
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the vortex doesn’t actually roll back toward the ground by reversing itself like a wheel, the 
density distributions show that the body forces actually push the primary vortex back toward the 
ground plane and the lateral density variations even push the vortex toward the symmetry plane. 
Clearly, density stratification can confine then destroy trailing line vortices much more rapidly 
than any other process considered in this study. 

Conclusions 

This study has shown that trailing line vortex flows can be resolved in the vicinity of a 
ground plane. Vortex trajectory comparisons with the experimental measurements of Liu and 
Smsky are in quite good agreement. The preliminary inclusion of modelled turbulence effects 
has not produced any startling results. On the other hand, density stratification can have a very 
pronounced effect on vortex trajectories. At the admittedly high stratification levels considered 
here, we have shown that vortex hazard alleviation could be facilitated by density stratification 
effects within the immediate vicinity of airport runways. 
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Appendix A. Asymptotic Solution for Flow Field Initialization 

Asymptotic expansions, in terms of e = y/t JW e , were used in both the outer flow and 
the inner flow. The two-dimensional, unsteady incompressible, Navier-Stokes equations were 
used as the governing equations in both regions. Employing van Dyke (1976) type matching 
procedures, the streamfunction and velocity, as well as pressure, were matched asymptotically. 
Analytic solutions were obtained when the inner solution was expanded further in time. 

In the following development, a two-term expansion is represented for this application. 
First, the expansion for streamfunction and velocity are presented. Then matching conditions 
are obtained by switching outer and inner variables according to van Dyke’s (1976) procedure. 
Since pressure can be derived if the velocity field is known in an incompressible flow, pressure 
matching is accomplished through the governing equations and will be developed after the 

governing equations in both regions are discussed. Finally, a closed form solution for the two- 
term expansions is developed. 

Expansions and Matching for Streamfunction and Velocity 

The appropriate streamfunction for the outer flow can be written 


* = *.(*.», <) + et,(x,!M) + 0[t J ], 


where 


and 


dV 0 d*i ^ rri 
dy dy +e dy + °[ C 1 


dV d * 0 dVx 21 

v --§i = -^-‘-dr+ 0 l‘ 2 l 


Similarly, the inner flow is represented as 


(A.1) 


(A.2) 


(A. 3) 


0 = 2 e[il> 0 (x,T),t) + e0i(x,»7,f)] + 0[e 3 ] 


(A.4) 



where 
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v = y/ 2c. 


Here, the inner velocity components are 


and 


50 

dy 


U = — = 2e 


1 50 o 


2 1 dy 
50q 50i 

dy + dy 


1 50i \ 

2 dy) 


+ 0[e 2 ] 


+ 0[e 2 j 


(A.5) 


V — 


dtp 

dx 




+ 0[ £ 3 ] 


(A.6) 


Matching Conditions 


Following van Dyke's (1976) matching theorem, we rewrite the outer expansion in inner 
variable while the inner expansion is rewritten in outer variable y. We then rewrite both the 
streamfunction and the velocity and match the corresponding temis. It is shown subsequently 

that using this matching procedure allows the streamfunction and velocity to be matched 
simultaneously. 


Rewriting the outer expansion in inner variables and expanding in terms of e , we have 


* = 2 e/ 7 , t) + €$!(*, 2ey, t) + 0[t 2 ] 


^o(x, 0, f) + t 


+ 0[e 2 ] 


(A.7) 


„ d<H 0 , 


2 V^rr-{x,0,t) + ^i(z,0,<) 


dy 2 


+ o[*l 


(A. 8) 


»0-e 


2 ^< x '°' , > + +o[e 2 ] 


(A.9) 
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Rewriting the inner expansion in outer variables and expanding in terms of e, 


iff = 2e 


+#1 ( x ’2?0j +0 f t3 l 


2e0 o (j, oo, t) + O [e 2 ] 


(A. 10) 


dib 0 / 

oo, f) + oo, f) + Oft 


(A.l 1) 


V % 


9 Mo ( 
dx ( ‘ r,0 °’ 


0 + O[e 2 ] 


(A. 12) 


From Eqns. (A. 10), (A. 11) and (A. 12), one can see that in the boundary layer the stream- 
function and y-component of velocity is one order lower in e than the x-component of velocity. 

Now we match streamfunction, x-component of velocity and y-component velocity, respec- 
tively, but to different order. We obtain 

0[c°] : 

*o(x, 0,0 = 0 (A. 13) 

— (x,°,t) = -^( X)00j i) = (A. 14) 

dty 0 

~d^(x,0,t) = Q (A.15) 


Of, 1 ] : 


d^ 0 

0 + 0 , t) 


= 2ip 0 (x,i h t) 


ff — * OO 


(A. 16) 


I f } — ►OO 


d 2 * 0 




2rj—( x ,0, t ) + -^ M t) 


Tj — OC 


dth 

djj 


(x,T},t) 


(A. 17) 


I TJ— >OC 
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Mo ( ., 


(A. 18) 


> *7—00 


As expected, Eqns. (A.13)-(A.18) either give the bottom boundary condition for the outer 
flow or the top boundary condition of the inner-flow. Equations (A. 13) and (A. 15) are equivalent 
and they just give one boundary condition for V 0 which is obviously the slip-free boundary 
condition. Equation (A. 14) is the top boundary condition for ip 0 which must match (at the edge) 
the x-component of velocity, U 0 . Equations (A. 16) and (A. 18) are also equivalent and they 
provide a boundary condition for while Eqn. (A. 17) is a boundary condition for The 
equivalence between Eqns. (A.13) and (A.15), (A. 16) and (A.18) shows that the streamfunction 
and velocity matching can be achieved simultaneously without contradiction in this procedure. 
The physical boundary conditions are thus retained by the mathematical manipulation. 

From the matching procedure, the solution for the two series can proceed sequentially; first 

the outer solution tf 0 , then the inner solution can be developed; followed by solution of 
^i, and so on. 


Governing Equations 

As stated previously, the governing equations used are the two-dimensional unsteady, incom- 
pressible Navier-Stokes equations. Since continuity is implied by introducing the streamfunction, 
we will concentrate on the momentum equations. 

Outer Flow 


Because viscous terms do not have any direct effect until terms of order e 2 are important, it 
is not necessary to include viscous effects for the two-term outer solution. Thus, we have 




dP 

dx 


(A. 19) 
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™ + v™. + 

dt dx dy dy 


(A.20) 


In addition, if we assume the flow is irrotational — especially near the edge of the boundary 
layer — the potential flow relation 


V 2 * = 0 


(A.21) 


can be used instead of Eqns. (A. 19) and (A.20). That has been done here and the Euler equations 
have been used only to determine the pressure field. 

For different orders of e, we get 


0[e°] : V 2 tf 0 = 0 


(A. 22) 


0{t x ] : V 2 tf 1= 0 


(A. 23) 


Since the outer solution here is only used for the purpose of generating the boundary-layer 
solution, the potential solution may or may not be applied to the whole outer flow field. 

Inner Flow 


The dimensionless component of the boundary layer momentum equation can be written 


du du du dp l d 2 u 

dt U dx V dy~~fa + fcfy 2 (A.24) 

when we have already assumed 0 and $£ can be neglected because both terms do not appear 
in the first two terms of the solution. From Eqns. (A.4) and (A.5), we get 


du 

~di 


gVo _ H d 2 0 o ^ 1 dtp i ( d 2 tpi 
dr)dt 2t df] 2 +e [Tt~d^ + d^dt 


V dVi 
2 1 dr i 2 


+ 0[e 2 ] 


(A.25) 



30 


and 


If we let 


du 

1 

’^ 2 V 0 


dy ' 

“ 27 

5t? 2 

+e <vJ 

d 2 u 

1 

■<9 3 !/) 0 

, d 3 0i" 

dy 2 = 

4e 2 

dr? 3 

v J 


+ 0[e] 


+ 0[t'] 


(A.26) 


(A. 27) 


P - Po(x, t) + epi(x, t) + 0[t 2 ], (A. 28) 

these expressions (Eqns. A.25-A.28) along with Eqn. (A.6) can be substituted into Eqn. (A.24), 
to obtain the equations for different orders of t: 

0 [<°] •' 

_ V Ho d 2 V>0 _ dhpodipo 

drjdt 2 1 drj 2 drj drjdx dr) 2 dx 

_ _dpo_ 1 dHo (A.29) 

dx 4 1 dr) 3 


The three boundary conditions are the non-slip boundary condition and matching conditions: 
(from Eqn. A. 14): 


i> 0 (x,o,t) = o 


(A.30) 


Ho 

dr) 


{x,0,t) = 0 


(A.31) 


Ho 

dr) 


= -~(x,0,t) = U o 


(A. 32) 


0 [* 1 ] : 

}_H\ Hh i)_dHi drp 0 d 2 ip l 

2 1 dr) dr/dt 2 1 dr) 2 + dy dydx 

+ Hi dHo _ HodHi dip i d 2 xp 0 

dr) drjdx dx dy 2 ~ ~dx dr) 2 

_ dpi + 1 jjVi 

da: 4f drj 3 


(A.33) 
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The three boundary conditions for xp { are no-slip conditions and matching condition 
Eqn. (A. 17): 


0l(x,O, t) = 0 


(A. 34) 


djn 

drj 


(x,0,0 = 0 


(A. 35) 


dip i 


T} — ♦ OO 




(A. 36) 


From Eqn. (A.22), we have 


d 2l $o d 

dx 2 °* ^ + ~dy^^ X ' °’ 0 = 0 


while Eqn. (A. 13) yields tf o (x,0,t) = 0. 
Hence 


a» a 1 ’ ’ ' a* a (x ’ () a ^ — 

Therefore, Eqn. (A. 36) can be rewritten: 


= 0 


(A.37) 


001 

drj 


= -^-(x,0,<) = f/i 


(A.38) 


Pressure Gradient Matching 

As mentioned earlier, to second order we have §£ = 0 in the boundary layer for the 
asymptotic expansions developed here. Consequently, we have 


dP _dp 
0 X y=o di 


(A. 39) 
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From the governing equation for the outer flow (Eqn. A. 19), we have 


dU 

T T dU 

dP 

dt 

+U q ~ 

y=0 OX 

y = 0 dx 


Since from Eqn. (A.2) 


(A.40) 


dU _ d 2 fr 0 d 2 $i l d'i 

dt ~ dydt + 6 d^dt + Yflhj ’ (A - 41 ) 


we obtain 

0[t°] : 


dUo dU 0 _ dp 0 
dt +Uo dx ~ dx 


(A.42) 


0[e ! ] : 


dU x 

dt 


+ h u ' + v ‘i £ + Ui 


dUo 

dx 


dpi 

dx 


(A.43) 


Solution for Inner Flow 

From Nam (1990), it is known that in an unsteady boundary layer, a two-series expansion 
is needed, in terms of e and t. In the following, we write each inner solution in different orders 
of e as a power senes in terms of t. Since we just use the solution as the initial condition in 
this problem, we can choose small t and then only the first several terms are needed. A closed 
form solution is thus obtained for the first two terms of the expansions. 

0[e°] : 

Let 

OO 

= X 0 

p = o 


(A.44) 



Then the first two terms are: 
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0o Uofoi +t(^ ^ fn + U 0 /i2^ + 0[t 2 ] 

(A.45) 

Substituting Eqn. (A.45) into Eqn. (A.29), with boundary conditions, Eqns. (A. 30), (A.31), 

(A.32) and employing the pressure gradient Eqn. (A.42), we obtain, for different orders of t: 

0{t° ] : 

C + 2 nfHx = o, 

subject to: 

(A.46) 

o 

II 

o 

II 

o' 

f-H 

(A. 47) 

/ij(oo) = 1 

(A.48) 

0[t'] : 

fu + 2 nfn ~ 4 fn = -4 + 4/J, 

subject to: 

(A. 49) 

O 

II 

o' 
1— 1 

II 

o' 

f-H 

(A.50) 

/ii(°°) = o 
and 

(A.51) 

/n + 2 7 /" 2 - 4/{ 2 = -4 + 4(/o,) 2 - 4/oi/oj 

subject to: 

(A.52) 

/h(0) = A ( 0)=0 

(A.53) 
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/j' 2 (oo) = 0 


(A. 54) 


The solutions for f 01 , f u and / 12 arc given in the text as Eqns. (1 1)-(13). 

0 [«‘] : 

Let 


Then the first two terms 


^1 = ^ tP 4>\{x,T],t) 

P= 0 


(A.55) 


4>i - + t(^j±g u + + u ^9u^j + 0[t 2 } (A.56) 

are gotten by substituting Eqn. (A.56) into Eqn. (A.33), with boundary conditions Eqns. (A.34), 

(A.35) and (A.36) and pressure gradient condition Eqn. (A.43). We obtain for different orders 
of t 

0[t°] : 

9oi + ~ 2<7oi = ~2 (A. 57) 

subject to 


0oi(O) = soi (0) = 0 


(A.58) 


0Ol(oc) = 1 


(A.59) 


0 [<‘] : 

9n + 2 n g'h - 6 g' u = -4 + 4^ 

subject to: 


(A.60) 


9u(0) = 0n(O) = 0 


(A.61) 
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Sn(oo) = 0 


and 


9\2 + 2 V 9 i 2 - & g'u = -4 + 4/oj^oj 


subject to: 


^ 12 ( 0 ) = g[ 2 (0) = 0 


5i 2 (oo) = 0 

and 

913 + 2 T) g' u - 6g u = -4 + 4 f' 01 g' Qi - 4 f 01 g'^ 

subject to: 


^ 13 ( 0 ) = ^ 3 ( 0 ) = 0 


9n(<x> ) = 0 


The solutions for^i, gu, g 12 and g 13 are given in the text as Eqns. (14)-(17). 


(A.62) 


(A. 63) 


(A. 64) 


(A.65) 


(A. 66) 


(A.67) 
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Appendix B. Derivation of Fluctuation Density Equation through Conservation of Mass 

In incompressible flow, the local density is a function of temperature and here it is assumed 
linearly related to temperature via 

P = Po[\ ~ a(T - T 0 )] (B .1) 

where a is the volumetric coefficient of thermal expansion and a << 1. 

We define dimensionless density and temperature as 

P=(p-Po)/p 0 (B.2) 

e = (T~ T 0 )/T 0 (B , 3) 

respectively. 

For the case of a linearly stratified ambient fluid, the quiescent ambient state, 0, is assumed 
given by 

9 = Py (B.4) 


Hence, the ambient density departure is 

P ~ ~ q (T — T 0 ) = —afiyT 0 


(B.5) 


and 


dp 

Ty = 


(B.6) 


Now, with the density fluctuation and temperature fluctuation defined as 


P = P-P 


(B.7) 
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O' = 6-9 

we have 


(B.8) 


P = —otT 0 9' 


(B.9) 


From the continuity equation for compressible flow 

dp dp dp 

di +U di + % +pV ' <: = 0 ’ (B.10) 

we obtain the dimensionless equation 

dp dp dp 

dt + U di + V ~fy + (1 + ^ )V ' “ = 0 (B.l 1) 

Substituting Eqn. (B.7) into Eqn. (B.l 1), we get 

dp 1 dp 1 dp' dp _ 

~dt +U 'dI + V 'd^ + V d^ + V ‘ il+ (p + p')v.u = 0 (B.12) 


From Eqns. (B.6) and (B.9), one can see that 

~ 0[a] 

ay 


p ~ 0[a] 

Hence, for terms of 0[a°], Eqn. (B.12) gives 

which recovers incompressibility. 

Then for 0[a], Eqn. (B.12) gives 

% + " + % = ~ v f v < B - 14 > 

which is the same as Eqn. (47) in the text, which was derived from the conservation of energy 


V • u = 0 


do' do' 


n 


(B.l 3) 


equation. 
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